Interatomic Methods for the Dispersion Energy Derived from the Adiabatic 
Connection Fluctuation-Dissipation Theorem 
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Interatomic pairwise methods are currently among the most popular and accurate ways to include 
dispersion energy in density functional theory (DFT) calculations. However, when applied to more 
than two atoms, these methods are still frequently perceived to be based on ad hoc assumptions, 
rather than a rigorous derivation from quantum mechanics. Starting from the adiabatic connection 
fluctuation-dissipation (ACFD) theorem, an exact expression for the electronic exchange-correlation 
energy, we demonstrate that the pairwise interatomic dispersion energy for an arbitrary collection of 
isotropic polarizable dipoles emerges from the second-order expansion of the ACFD formula. More- 
over, for a system of quantum harmonic oscillators coupled through a dipole-dipole potential, we 
prove the equivalence between the full interaction energy obtained from the Hamiltonian diagonal- 
ization and the ACFD correlation energy in the random-phase approximation. This property makes 
the Hamiltonian diagonalization an efficient method for the calculation of the many-body dispersion 
energy. In addition, we show that the switching function used to damp the dispersion interaction 
at short distances arises from a short-range screened Coulomb potential, whose role is to account 
for the spatial spread of the individual atomic dipole moments. By using the ACFD formula we 
gain a deeper understanding of the approximations made in the interatomic pairwise approaches, 
providing a powerful formalism for further development of accurate and efficient methods for the 
calculation of the dispersion energy. 



I. INTRODUCTION 

Van der Waals (vdW) forces are ubiquitous in nature, 
and they play a major role in determining the struc- 
ture, stability, and function for a wide variety of systems, 
including proteins, nanostructured materials, as well as 
molecular solids and liquids. A significant attractive part 
of the vdW energy corresponds to the dispersion energy, 
which arises from correlated fluctuations between elec- 
trons. Therefore, accurate treatment of the dispersion 
energy is essential for improving our understanding of 
biological and chemical systems, as well as (hard and 
soft) condensed matter systems in general. Many en- 
couraging ideas and methods have been proposed in re- 
cent years for approximately including the missing long- 
range dispersion interactions in density functional theory 
(DFT) pHl. Despite significant progress in the field of 
modeling vdW interactions during the last decade, many 
questions still remain unanswered and further develop- 
ment is required before a truly universally applicable 
(accurate and efficient) method emerges. For example, 
interatomic vdW potentials are frequently employed for 
the modeling of hybrid inorganic/organic interfaces @- 
Hoj . neglecting the relatively strong long-range Coulomb 
screening present within inorganic bulk materials. On 
the other hand, the popular non-local vdW-DF function- 
als [ITl - fl3j use a homogeneous dielectric approximation 
for the polarizability, which is not expected to be accu- 
rate for molecules. Nevertheless, the interaction energies 
between small organic molecules turn out to be reason- 
ably accurate. Understanding the physical reasons of 
why these different approaches "work well" outside of 
their expected domain of applicability is important for 
the development of more robust approximations. 



Interatomic pairwise dispersion approaches, typically 
dubbed DFT-D Q or DFT+vdW 0, are among the 
most widely used methods for including dispersion energy 
in DFT. Such approaches approximate the dispersion en- 
ergy in a pairwise fashion, i.e. as a sum over unique 
atom pairs. Despite their simplicity, these effective pair- 
wise models provide remarkable accuracy when applied 
to small molecular systems, in particular when accurate 
dispersion coefficients (C n ) are employed for atoms in 
molecules 15]. Only recently have efforts been focused 
on going beyond the effective pairwise treatment of vdW 
contributions, for example, the role of the non-additive 
three-body interatomic Axilrod-Teller-Muto term was as- 
sessed (Ha [III- Furthermore, an efficient and accurate 
interatomic many-body dispersion (MBD) approach to 
dispersion interactions has recently been proposed [l8j . 
The MBD description of vdW interactions is essential 
for the description of extended molecules and molecular 
solids; however, the influence of MBD interactions can 
already become significant when considering the interac- 
tions between relatively small organic molecules [ill [lj| . 

Despite the popularity and the relative accuracy of the 
DFT-D and DFT+vdW methods, they are still widely 
perceived to be based on ad hoc assumptions. For the 
dispersion interaction between two spherical atoms i and 
j, the pairwise C$ RZ formula has been known since the 
seminal work of London [20l[2l|. However, to the best of 
our knowledge, the generalization of London's formula for 
an arbitrary collection of N spherical atoms has not been 
explicitly derived from first principles. Furthermore, at 
short interatomic distances, the dispersion interaction is 
significantly weaker than the corresponding asymptotic 
expansion, and ad hoc approximations become necessary 
for the functional form of the damping. It is shown in this 
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work that the switching function used to damp the dis- 
persion interaction at short distances arises from a short- 
range screened Coulomb potential, whose physical role is 
to account for the spatial spread of the individual atomic 
polarizabilities. 

Since the dispersion energy arises from correlated fluc- 
tuations between electrons, it is intrinsically a many- 
body phenomenon, an accurate description of which re- 
quires a quantum mechanical treatment. For this pur- 
pose, the adiabatic connection fluctuation-dissipation 
(ACFD) theorem, which provides a general and exact 
expression for the exchange-correlation energy [22|, [23| . 
allows us to calculate the dispersion energy in a seamless 
and accurate fashion, naturally including many-body ef- 
fects. In this work, we show that the ACFD theorem pro- 
vides a firm theoretical basis for the development and un- 
derstanding of interatomic pairwise and many-body dis- 
persion methods. In particular, we derive the well-known 
Ce/R 6 interatomic pairwise summation formula as the 
second-order expansion of the ACFD correlation energy, 
and demonstrate that this formula is valid for an arbi- 
trary collection of N fluctuating dipoles, each of which 
is characterized by an individual frequency-dependent 
polarizability. By applying the ACFD formalism we 
also prove, for a system of quantum harmonic oscillators 
(QHOs) coupled within the dipole approximation, the 
mathematical equivalence between the exact dispersion 
energy and the correlation energy in the random-phase 
approximation (RPA). This analytical result makes the 
coupled-oscillator model (with a relatively minimal com- 
putational cost) an ideal candidate for the inclusion of 
MBD effects in DFT. Finally, we show the relevance of 
MBD energy on the binding energies of dimers in the S22 
database. The full many-body description consistently 
reduces mean relative errors with respect to the inter- 
atomic pairwise approximation, showing the largest im- 
provements for the most extended systems. The ACFD 
formula leads to a deeper understanding of the approx- 
imations made in the development of the DFT-D and 
DFT+vdW approaches, and provides a natural formal- 
ism for further improvement of methods for computing 
the dispersion energy. 



II. THE PAIRWISE INTERATOMIC 
DISPERSION ENERGY FROM THE ACFD 
FORMULA 

The ACFD theorem is an exact expression for the 
exchange-correlation energy of a system of nuclei and 
electrons, described by a response function x(r, r', u>) [22I 
[23|]. The response function x measures the response at 
point r due to a change of the potential at point r' as 
a function of time or (Fourier-transformed) frequency u. 
Here our interest lies in the dispersion energy, which is 
contained in the electron correlation energy. Therefore 
we start by writing the ACFD formula exclusively for 
the correlation energy (Hartree atomic units are used 



throughout) 



2tt 



duj / d\Ti[( X x-Xo)v], (1) 



in which xo is the bare or non-interacting response func- 
tion, which can be computed given a set of single-particle 
orbitals [24, 25] fa (with corresponding eigenvalues and 
occupation numbers /j) as 

X0 {r,r',u) =2 2JJi-fi) J , (2) 

X\ is the interacting response function at Coulomb cou- 
pling strength A, defined via the self-consistent Dyson 
screening equation x\ = Xo + Xo(Av + f£ c )X\, « = 
|r — r'| _1 is the Coulomb potential, and Tr denotes the 
trace operator over spatial variables r and r'. Using the 
ACFD formula, the adiabatic connection between the 
non-interacting system (with A = 0) and the fully in- 
teracting system (with A = 1), yields the full correlation 
energy of the system of interest. Obviously, the corre- 
lation energy obtained from the ACFD formula contains 
the full many-body dispersion energy as well as other 
electron correlation effects. 

In practice, the exact form of the exchange-correlation 
kernel fg c in the Dyson equation is not known. Neglect- 
ing the explicit dependence of f xc on A, analytic inte- 
gration can be carried out over A. This is the case, for 
example, for the widely employed random-phase approx- 
imation (RPA) [26l . |27( , or the full potential approxima- 
tion (FPA) In the RPA, f xc = 0, while in the FPA 
Xx = Xii *- e - the A integration is carried out using the 
fully interacting response function. In what follows, we 
will employ the RPA method, which has been shown to 
yield reliable results for a wide variety of molecules and 
extended systems [29U43I ] . Using the Dyson equation, the 
ACFD correlation energy expression in Eq. ([T]) takes on 
the following form in the RPA 
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1 - Xvxo 



(3) 



Integration over A in Eq. Q leads to the following ex- 
pansion for the correlation energy in terms of xo v |44.|45| 
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n=2 



-Tr[(xo«) 



(4) 



Let us now apply the ACFD- RPA approach to a collec- 
tion of fluctuating dipoles representing the atomic system 
of interest. Each atom i is characterized by its position 
fi = {xi,yi,Zi} and a frequency-dependent dipole polar- 
izability oei(kj). For the moment, we assume that the 
atoms are separated by a sufficiently large distance, al- 
lowing us to use the bare Coulomb potential to describe 
the interaction between dipoles. The general case will be 
addressed in the next section. The response function for 
each atom i takes the form [46| 

Xi(r, r', iu) = -a,(iu;)V r £ 3 (r - r ; ) ® V r ^ 3 (r' - r ; ), (5) 
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where 5 is the three-dimensional Dirac delta function, 
and (g> is the tensor (outer) product. 

Now it will be shown that the second-order (n = 2) 
term in Eq. Q yields the pairwise interatomic dispersion 
energy. In the remainder of this section we drop the A 
index of the response function, because the atomic polar- 
izabilities can be derived from a mean-field or an explicit 
many-body calculation. Furthermore, the second-order 
n = 2 term in Eq. (Q| turns out to be the same if the 
expansion is carried out in terms of xo or xi l43|- For a 
collection of N atoms in the dipole approximation, the 
XV matrix can be written as AT. Here, A is a diagonal 
3iV x 3N matrix, with —cti(iui) values on the diagonal 
blocks. The T matrix is the dipole-dipole interaction 
matrix, with 3x3 ij tensors given by Ty = V ri ® V rj v%j 
(Ta — 0). For two atoms (ai(iw) and 0:2 (iw)) separated 
by a distance R = \r\ — ra| on the x axis, the AT matrix 
becomes 
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(6) 



The second-order (n = 2) term of the ACFD-RPA ex- 
pression in Eq. (j4)) with the above matrix as input leads 
to 

1 f°° C 12 
E c%pa = -^ / dwai(iw)a 2 (iw)Tr[Ti2 2 ] = 

(7) 

where the Casimir-Polder identity has been used to de- 
termine Cg 2 from the corresponding dipole polarizabil- 
ities. The above equation is the familiar expression for 
the long-range dispersion interaction between two atoms. 

Equation [7] can be easily generalized to the case of N 
atoms. The use of the trace operator in Eq. Q requires 
multiplication of column i of the AT matrix by the corre- 
sponding row i. Since for any given TJy, Tr[Ty 2 ] = 6/Rfj, 
where Rij is the distance between atoms i and j, the 
second-order expansion of Eq. (|4]) for N atoms leads to 
this simplified form 

i j y 

This of course is the familiar expression for the dispersion 
energy for an assembly of N atoms used in the DFT-D 
and DFT+vdW methods. We note that the above deriva- 
tion of the pairwise dispersion energy from the ACFD 
formula does not make any assumptions regarding the ge- 
ometry of the atomic assembly or the functional form of 
the frequency-dependent polarizabilities. Furthermore, 
we note that employing the FPA instead of the RPA 
would not change the conclusions presented above. 



While the second-order expansion of Eq. (|4|) for 
isotropic polarizabilities yields the familiar pairwise in- 
teratomic dispersion energy given by Eq. (j8|), the former 
equation is more general. It allows for the use of full po- 
larizability tensors, enabling an anisotropic treatment of 
the dispersion energy. In this regard, the polarizability 
anisotropy has been found to play a non -neg ligible role 
for intermolecular dispersion interactions [471 . |48( . 

We note that the higher-order terms in the RPA ex- 
pansion of the correlation energy include two contribu- 
tions: higher-than-pairwise many-body interactions (up 
to A-th order) and the response screening (up to infinite 
order). As an example of the beyond-pairwise many- 
body interactions captured in the RPA expansion of 
the correlation energy, the third-order term includes the 
well-known Axilrod- Teller- Muto three-body energy [49j . 
The higher-order response screening can be easily illus- 
trated for two interacting atoms i and j by expanding 
Eq. (@} (explicit dependence of the polarizabilities on iui 
assumed) : 

1 f°° , / 6a,a 7 36afa 2 \ 

e ^ a= -tJ q d "[-w + ^w L ---j- (9) 

in which the second-order term corresponds to the "stan- 
dard" Ce/R e dispersion interaction between i and j, and 
the higher-order terms (which only survive with even 
powers of n) correspond to the response screening of the 
polarizability of atom i by the presence of atom j and 
vice versa. Further analysis of these many-body contri- 
butions is presented in Sections IV and V. 

III. THE DAMPING OF THE DISPERSION 
ENERGY AT SHORT DISTANCES 

Correlation energy calculations carried out using the 
ACFD formula usually employ the response function xo 
computed using a set of occupied and virtual one-particle 
orbitals [see Eq. @], determined from semilocal DFT, 
Hartree-Fock, or hybrid self-consistent-field calculations. 
In this scenario, xo is typically a fairly delocalized object, 
which includes orbital overlap effects between occupied 
and virtual states. However, even in this case, the use of 
certain approximations for the exchange-correlation ker- 
nel fxc can lead to divergencies for close inter-particle 
separations [30] . When the response function is localized, 
leading to a diagonal form, the details of the overlap be- 
tween orbitals are lost. For example, this is clearly the 
case for an assembly of fluctuating point dipoles. When 
two point dipoles come into close contact, the Coulomb 
interaction between them diverges. In fact, depending 
on the absolute values of the polarizability, the head-to- 
tail alignment between two dipoles can lead to an infinite 
polarizability even for a finite (non-zero) separation be- 
tween the dipoles [5(| ■ This is clearly an unphysical situ- 
ation, which is mitigated in practice by the finite extent 
of electronic orbitals. From a slightly different perspec- 
tive, the dipole moment should be spread out in space, 
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and this effect naturally introduces damping when the 
polarizabilitics overlap. 

The most straightforward way to incorporate overlap 
effects for a set of fluctuating point dipoles is through 
a modification of the interaction potential. Thus, in- 
stead of using the bare Coulomb potential v = |r — r'| _1 , 
a modified potential should be employed that accounts 
for orbital overlap at short distances. We take the po- 
larizability of each atom i to correspond to a quantum 
harmonic oscillator (QHO) 



1 + (w/w ) 2 



(10) 



where ao is the static dipole polarizability, and luq is 
an effective excitation frequency. Since the ground state 
QHO wavefunction has a gaussian form, the interaction 
between two QHOs (or atoms) i and j leads to a modified 
Coulomb potential 



ed(rij/aij) 



where r« is the interatomic distance, cry is an effective 
width, (Jij = \/o~f + o~? obtained from cr, and <7 7 , the 



gaussian widths of atoms i and j , respectively. Since the 
polarizability relates the response of a dipole moment to 
an applied electric field, the <7j and Uj parameters cor- 
respond to the distribution of the dipole moment, and 
not of the charge. The width of the gaussian distribu- 
tion is directly related to the polarizability in classical 
electrostatics |5lj . 

Using Eq. (JTTJ) , the dipole interaction tensor for atoms 
i and j becomes 



rpab 



3r a r 6 - r?-(* n6 ( ) 



saturates to unity, typically for distances 20% larger than 
the sum of the van der Waals radii of the two atoms. 
Besides these two constraints the functional form of the 
damping is essentially arbitrary. Some evidence sug- 
gests [53[ that the binding energies are not significantly 
affected by the functional form of the damping, as long as 
two adjustable parameters are used. One of these param- 
eters controls the steepness, while the other determines 
the onset of the damping in terms of the distance. Typi- 
cally, only one of these parameters is adjusted for a given 
DFT functional by minimizing the error with respect to 
high-level quantum-chemical binding energies. Another 
disadvantage of the damping function, as used in DFT-D 
and DFT+vdW methods, is the need to define van der 
Waals radii, which are not quantum mechanical observ- 
ables. 

Inspection of Eq. (TT2"j) shows that, when using a QHO 
approximation for the spatial spread of the polarizabil- 
ity, the damping function is more complicated than a 
purely multiplicative function. This complication arises 
(11) due to the last exponentially decaying term in . Al- 
though this conclusion is based on a QHO model, the 
same conclusion holds for other models, such as hydro- 
genic atoms. We conclude that our findings are likely 
to be valid in general, meaning that the damping func- 
tion must be derived from a model Coulomb potential 
that naturally accounts for short-range dipole distribu- 
tion overlap effects [18j . The coupling of the dispersion 
energy to a given DFT functional might still require em- 
pirical parameter(s), as we illustrate below. However, a 
seamless coupling may be achieved by using the range- 
separation of the Coulomb potential in the calculation of 
the DFT correlation energy. The approach presented in 
this section can be employed in the development of such 
a range-separation procedure. 



y 

1 r a r b 



3 2 

ij ij 



(12) 



where a and b specify general Cartesian coordinates 
(x,y,z), r a and r b are the respective components of the 
interatomic distance ry, and <$y is the Kronecker delta 
function. It can be clearly seen that the above expression 
reduces the interaction between dipoles at short distance, 
in comparison to the bare dipole interaction potential. 
Even in the zero-distance limit, it converges to a finite 
value. Therefore, the description of the polarizabilitics 
by a dipole distribution instead of a point naturally in- 
troduces short-range damping effects, which have been 
so far included using ad hoc models in the DFT-D and 
DFT+vdW approaches. 

Both DFT-D and DFT+vdW methods use a distance- 
dependent damping function /d a mp( r y')> which multiplies 
the Cq R~j 6 dispersion energy. The function /damply) 
converges to zero or a small finite value at zero distance 
between two atoms [52| . At large distances, /damply) 



IV. EFFICIENT EVALUATION OF THE 
ACFD-RPA ENERGY FOR A SYSTEM OF 
QUANTUM HARMONIC OSCILLATORS 

The ACFD-RPA approach to the correlation energy 
has proven to be promi sing for a wide variety of molecules 
and extended systems [29l443| . The largest drawback of 
ACFD-RPA calculations is their relatively high compu- 
tational cost, resulting in a steep increase in the required 
computational time with system size. The conventional 
scaling of the ACFD-RPA calculations is N 5 , where N 
is the number of basis functions. This can be reduced 
to N , when resolution-of-the- identity, or density- fitting, 
techniques are employed to compute the four-centered 
two-electron Coulomb integrals [45j . 

If we only aim at computing the long-range vdW en- 
ergy, it is possible to associate a single quantum harmonic 
oscillator (QHO) to every atom. For such a system, in 
which the QHOs interact within the dipole approxima- 
tion, one can circumvent evaluation of the four-centered 
two-electron Coulomb integrals (and costly summations 
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over the virtual/unoccupied functions) via application 
of the QHO selection rules [54T - I561 ] . In this section, we 
will demonstrate that the ACFD-RPA correlation energy 
for a system of QHOs interacting through a dipole po- 
tential is equivalent to the interaction energy obtained 
from diagonalizing the corresponding Hamiltonian ma- 
trix. Within this formalism, the computational scaling 
will be reduced to iV 3 (via the diagonalization step), 
where N is simply the number of atoms in the molec- 
ular system of interest. Further computational savings 
can be obtained if one computes the interaction energy 
with efficient path integral techniques [13, [Hj], instead 
of diagonalizing the coupled QHO Hamiltonian matrix, 
allowing for the efficient computation of the many-body 
dispersion energy during Monte Carlo (MC) and molec- 
ular dynamics (MD) simulations. 

Throughout the remainder of this Section, we will 
restrict the derivation to a system of one-dimensional 
QHOs in order to simplify the notation. The extension 
to three-dimensional QHOs is straightforward and does 
not alter the conclusions. 



A. The Hamiltonian for a System of Coupled 
QHOs 

For a system of iV QHOs interacting within the dipole 
approximation, the Hamiltonian can be written as [181 ]: 



H = - 



^V 2 



Hp 



p 



N 2f2 

2 

p 



N 



p>q 

where £ p represents the displacement of the p-th oscilla- 
tor from its equilibrium distance weighted by the square 
root of its mass. The first two terms in this Hamiltonian 
represent the kinetic energy and the confining potential 
corresponding to a set of independent QHOs with unit 
charge and characteristic frequency lu p (i.e., the eigenval- 
ues of the non-interacting, single-particle QHO Hamilto- 
nian matrix). The third term describes the interoscilla- 
tor coupling via the dipole-dipole interaction, which also 
depends on lu p , as well as a p , the QHO static dipole po- 
larizability, and T pq , the N x N dipole-dipole interaction 
tensor (as defined in Section II). Due to the quadratic 
(bilinear) dependence of the Hamiltonian on the £ p co- 
ordinates, it is possible to obtain an exact solution for 
the QHO interaction energy via diagonalization of the 
following cQ HO matrix: 

Cpf 10 = S pq^l + C 1 - Spq^uipUg^/a^Tpq . (14) 

The resulting interaction energy, -E Cj qhOj is then com- 
puted as the difference between the eigenvalues of the 
coupled system of QHOs (obtained via diagonalization of 
the C QHO matrix) and the eigenvalues of the uncoupled 
system of QHOs (the characteristic frequencies), i.e., 



E, 



e,QHO 



2 ^ 



From Eq. (|15[) . the interaction energy for a system of 
coupled QHOs is characterized by a set of normal modes, 
whose corresponding frequencies are shifted with respect 
to the non-interacting characteristic frequencies due to 
the presence of dipole-dipole coupling: 



(16) 



B. 



The ACFD-RPA Correlation Energy for a 
System of Coupled QHOs 



After performing the integration over the coupling con- 
stant A in Eq. ([3]), the ACFD-RPA correlation energy 
can be written in the following form (an alternative yet 
equivalent expression to Eq. (g|) [591 ]: 



RPA 



f duTr\hx(l-Xov)+Xov]. (17) 
^ Jo 



As discussed in Section II, xo v corresponding to a set 
of coupled QHOs can be written in matrix form as AT 
(See Eq. ©) by utilizing the QHO selection rules and 
the inherent locality of the QHO polarizabilities. Since 
Tr[AT] = 0, the ACFD-RPA correlation energy for a 
system of coupled QHOs can be written as 



£ c ,rpa = — / du: In det[C RPA 
27T Jo 



M] > (18) 



in which the C RPA matrix is defined as follows: 
C* PA M = 5 pq + (1 - 6 pq )a p (iw)T pq . 



(19) 



C. Deconvolution of the lo v and uj p Dependencies in 
the C RPA Matrix 

In order to disentangle the u>i and uii dependencies in 
the (7 RPA matrix, we seek to rewrite (7 RPA as the product 
of two diagonal w-dependent matrices which separately 
contain uj p and ui p , respectively. This is accomplished by 
first extracting the free QHO polarizabilities in the (7 RPA 
matrix via 



C RPA (k;) = -A(iu)B(iu) 



(20) 



where A pq {iui) = — 6 pq a q (iLu) (c.f. Section II) is a diago- 
nal N x N matrix which only depends on the uncoupled 
characteristic frequencies, u) p , and 



B pq (\Lo) = S pq a p (iuj) 1 + (1 - d~ P q)T pq 



(21) 



is a non-diagonal N x N matrix which depends on both 
Lj p and Hi p . Since the dipole-dipole interaction tensor 
is frequency independent, one can follow the procedure 
suggested in Ref. and recast the B{\ui) matrix as 



( 15 ) B(iu) = B(0) + Dui 1 = B(0) + <5 OT (a p (0)a$-V 



(22) 



6 



in which the D matrix has been defined explicitly. In 
this form, B{ilo) can now be written in terms of a diago- 
nal matrix by solving the following generalized eigenvalue 
problem: 

B(0)t n = ui 2 n Dt n , (23) 

where t„ (n — 1,...,N) is a complete set of TV- 
dimensional vectors which diagonalize the D matrix, 
i.e., tf n Dt n — S mn . This is easily seen if one defines 
T = [ti,ta, • • • ,tjv], from which 

f T B{iuj)f = fi(iw) = S pq (uj 2 + Cj 2 p ) . (24) 

The (7 RPA matrix can now be written in terms of the two 
diagonal matrices A(iw) and f2(iw): 

CrpaM = -A(iu)(f T )~ 1 n(itu)f - 1 . (25) 

The matrix A(iw), containing the free polarizabilitics, 
will only depend on the uncoupled QHO frequencies u p . 
The coupled QHO frequencies, ui p , on the other hand, 
will be present inside the f2(iw) matrix through its de- 
pendence on the uj p . The relationship between the oj p 
and uj p can be determined by defining a diagonal matrix 
F such that F T F = D, 

F m = MM0K)-V2 (26) 

and observing the fact that B(0) = F T C QHO F, which 
allows for diagonalization of the coupled C^ HO matrix, 
i.e., 

f T B{0)f = fT F T c QKO F f _ q(q^ _ 

Hence, 17(0) is the matrix of the coupled eigenvalues of 
CQ HO and Q? = u>l 

D. Frequency Integration and the Contributions 
from the oj p and lo v Poles 

Before proceeding to the integration over frequency in 
Eq. (|18p. we first consider the logarithm operating on 

£<RPA. 

lndet[C RPA ] = hxdet[-An] + In det[(f T )~ 1 f ^} , (28) 
in which Eq. ([25]) was used. Since {ff T )~ 1 = D and 

(-A(iu;)Q(iu;)) M = 6 pq i (G$ + ^) , (29) 

all that remains in the expression for E c ^pa is N inte- 
grals of the form 

i ^ r frt + tjA , 



which, after integration by parts and the use of Eq. (|16l) , 
becomes 




This integrand shows both a pair of coupled (cj = ±iCu p ) 
and uncoupled (u> — ±iuj p ) QHO poles. Extending the 
integral to — oo by symmetry and closing the integration 
path in the upper imaginary half plane results in the fact 
that only the QHO poles possessing a positive imaginary 
component provide a non-zero contribution to -E c ,rpa- 
By explicitly performing the frequency integration, the 
coupled and uncoupled poles will contribute with a |wi|/2 
and —uii/2 term, respectively. From Eq. (|15p . the sum 
of these contributions for each of the TV QHOs yields 

Hence, the ACFD-RPA correlation energy for a set of 
QHOs coupled through a dipole-dipole potential is equiv- 
alent to the interaction energy that one obtains upon ex- 
act diagonalization of the Hamiltonian for a system of 
coupled QHOs. The coupled QHO normal modes nat- 
urally include many-body effects in complete analogy to 
the ACFDT-RPA energy. Although Eq. |15]) gives no fur- 
ther insight concerning the nature of these many-body ef- 
fects, we have demonstrated how the coupled QHO model 
naturally provides beyond-pairwise many-body energy 
contributions (up to TV-th order) and the RPA response 
screening (up to infinite order). As a result, diagonal- 
ization of the coupled QHO Hamiltonian allows for an 
effective RPA treatment of long-range vdW interactions 
at a significantly reduced computational cost. As such, 
the coupled QHO model represents a highly efficient and 
tractable method for the calculation of the many-body 
vdW energy in large scale systems. Furthermore, we 
stress that the present results do not depend on the choice 
of the T matrix. Any shape of the dipole-dipole interac- 
tion (as long as it remains frequency-independent) would 
not alter the validity of these conclusions. 



V. THE IMPORTANCE OF MANY-BODY 
EFFECTS FOR THE INTERATOMIC 
DISPERSION ENERGY 

The pairwise dispersion energy is an approximate form 
that we derived starting from the exact ACFD formula 
for the correlation energy in Eq. (TT|) . This derivation was 
based on two approximations: (i) Analytic integration 
over the adiabatic connection parameter A using f xc — 
(RPA) (ii) Second-order truncation of the logarithmic se- 
ries expansion resulting from (i). While the first approx- 
imation was proven to hold exactly for a system of cou- 
pled QHOs, the logarithmic series truncation limits the 
validity of the pairwise approximation to second-order 
perturbation theory. In this section, we assess the effect 
of going beyond the pairwise approximation. This anal- 
ysis is complementary to our recent work showing the 
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importance of many-bod y d ispersion energy for a variety 
of molecules and solids [l^, The difference here is 
that our analysis is now based on the rigorous ACFD- 
RPA expression. 

We model each atom i as a QHO as explained in the 
previous section. The input parameters ao and ujq in 
Eq. (I10p are obtained from first-principles by using the 
Tkatchenko-SchefHer (TS) method [l4[. In practice, for a 
given molecule, we carry out a DFT calculation using the 
PBE exchange-correlation functional 61]. The resulting 
self-consistent electron density is then partitioned into 
atomic contributions using the Hirshfeld approach f62| . 
Both ao and ujq are then defined as functionals of the 
atomic electron density. We note in passing that using 
densities from different functionals leads to essentially 
the same final results (for further details see Ref. [H]). 
The width Oi of every QHO in Eq. (fT0"| is determined by 

using Eq. (|12p in the limit of zero-distance between two 

i 

dipoles, <Ji 

For small and medium size molecules, the TS method 
yields intermolecular Cq coefficients in excellent agree- 
ment with experimental values obtained from dipole- 
oscillator strength distributions [14J. In this case, the 
input TS polarizabilities a correspond more closely to 
the fully interacting response function xi than to the 
non- interacting response function xo used in Eq. (U). Ex- 
panding ACFD-RPA formula in terms of the interacting 
response function xit instead of xo, leads to the following 
expression for the correlation energy |44j 

E c = -i- J™ du f>l)" (l - i) Tr[( X ivn (32) 

Comparing Eq. (|32|) to Eq. ((4]), it is clear that the 
second-order term remains the same, albeit operating 
now on xi instead of xo- This fact demonstrates that 
the origin of the polarizabilities does not modify the pair- 
wise additive expression for the dispersion energy, further 
strengthening our derivation of the pairwise interatomic 
methods from the ACFD formula in Section II. In con- 
trast, the coefficients of the higher-order terms differ sig- 
nificantly between Eq. (|32|) and Eq. (j4}. In fact, the 
terms with even powers of n carry a negative sign in 
Eq. (|3"2"j) . while all the terms in Eq. (Q} are positive. We 
note that the modified Coulomb potential W proposed 
in Ref. [l8[ leads to a similar switch in the sign of the 
many-body energy contributions. 

Widely used (semiVlocal and hybrid functionals in 
DFT, such as PBE [gy, PBEO [gi,^, and B3LYP [67], 
HH, are relatively successful for the short-range correla- 
tion energy. In contrast, our approach based on the QHO 
model, is constructed to accurately describe the long- 
range correlation energy. A seamless connection between 
a given DFT functional and the QHO model requires 
an explicit modification of the DFT functional correla- 
tion hole. This offers an interesting direction for future 
work. Here, instead we introduce an empirical parameter 



TABLE I. Performance of different functionals with disper- 
sion energy on the S22 database of intermolecular interac- 
tions. The errors are measured with respect to basis-set 
extrapolated CCSD(T) calculations of Takatani et al. [63| ]. 
Mean absolute relative errors (MARE in %), standard devi- 
ation (SD in kcal/mol), and mean absolute errors (MAE in 
kcal/mol) are reported. The postfix "-2D" means that the 
dispersion energy is added using the second-order expansion 
of Eq. (|32|l . while "-00D" means that the dispersion energy is 
computed to infinite order. The DFT calculations have been 
carried out with the FHI-aims code [64|] using a large numeric 
tier 3 basis set. For DFT, this basis set is converged to better 
than 0.05 kcal/mol compared to the basis set limit [TEI ] . 







MARE 


SD 


MAE 


PBE-2D 


2.06 


11.6% 


0.84 


0.64 


PBE-00D 


2.50 


7.1% 


0.57 


0.43 


PBE0-2D 


2.12 


11.4% 


0.96 


0.72 


PBEO-00D 


2.52 


7.2% 


0.66 


0.52 



TABLE II. Performance of PBEO functional including effec- 
tive pairwise dispersion (PBE0-2D) and the full many-body 
dispersion (PBEO-coD) on the dispersion-bound complexes 
contained in the S22 database, with respect to the basis-set 
extrapolated CCSD(T) calculations of Takatani et al. [gij. 
All values are reported in kcal/mol. The DFT calculations 
have been carried out with the FHI-aims code [6J] using a 
large numeric tier 3 basis set. For DFT, this basis set is con- 
verged to better than 0.05 kcal/mol compared to the basis set 
limit [TJ]. 



PBE0-2D PBEO-00D CCSD(T) 



Methane dimer 


-0.63 


-0.58 


-0.53 


Ethene dimer 


-1.64 


-1.37 


-1.48 


Benzene-Methane 


-1.45 


-1.48 


-1.45 


Benzene dimer (C211) 


-1.80 


-2.50 


-2.62 


Pyrazine dimer 


-3.00 


-3.23 


-4.20 


Uracil dimer 


-8.80 


-9.57 


-9.74 


Indole-Benzene Stack 


-3.15 


-4.54 


-4.59 


Adenine-Thymine Stack 


-9.86 


-11.80 


-11.66 



ft that multiplies the QHO-QHO interaction parameter 
a in Eq. (jTTJ) . A value of j3 larger than unity corresponds 
to an interaction that is shifted to larger distances, effec- 
tively capturing only the long-range part of the correla- 
tion energy. 

In order to assess the accuracy of different approxima- 
tions to the ACFD formula, we have chosen to use the 
S22 database of intermolecular interactions [69j . a widely 
used benchmark database with binding energies calcu- 
lated by a number of different grou ps using high-level 
quantum chemical methods [63l l69l|. In particular, we 
use the recent basis-set extrapolated CCSD(T) binding 
energies calculated by Takatani et al. (63[. These bind- 
ing energies are presumed to have an accuracy of ~ 0.1 
kcal/mol (1% relative error). This level of accuracy al- 
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lows an unbiased assessment of approximate approaches 
for treating dispersion interactions. Table U summa- 
rizes the results of our calculations on the S22 database. 
We used two different non-empirical DFT functionals: 
PBE [63 and PBEO [gHH. For every functional, we 
have carried out two tests: (i) using the second-order ex- 
pansion of Eq. (|52"j) . and (ii) using the full (infinite) series 
in Eq. (I32p . The /3 parameter has been adjusted for every 
combination of functional and method. We note that the 
approach presented here does not require the definition 
of van der Waals radii — all the necessary information is 
contained in the input frequency-dependent polarizabili- 
ties and the adjusted /3 parameter. 

An analysis of the performance of different methods on 
the S22 database in Table U reveals that the addition of 
the pairwise dispersion energy in the QHO approxima- 
tion yields a mean absolute error (MAE) with respect to 
CCSD(T), which is typically a factor of 2 larger than for 
DFT-D and DFT+vdW methods [lEGJ]- We note that 
the MAE can be reduced significantly by using a steeper 
functional form for the damping of the short-range in- 
teraction [l8j . The QHO approximation describes ev- 
ery atom as a single gaussian function, which leads to a 
smooth damping of the dispersion energy at short dis- 
tances. Thus, even hydrogen-bonded systems are sig- 
nificantly stabilized. Interestingly, the inclusion of the 
infinite order dispersion energy beyond the pairwise ap- 
proximation noticeably reduces the errors and increases 
the /3 value for every tested functional. Both of these 
results are desirable, since a larger value of j3 means that 
the dispersion energy is shifted to larger distances. In 
particular, the error is reduced for all dispersion-bound 
systems when going from PBE0-2D to PBEO-ooD as 
shown in Table HTl (the exception is the benzene-methane 
dimer, where both PBE0-2D to PBEO-ooD yield essen- 
tially the same results). The most pronounced devia- 
tion between PBE0-2D and PBEO-ooD is observed for 
the largest dispersion-bound adenine-thymine complex. 
The CCSD(T) method yields a binding energy of -11.7 
kcal/mol, while PBE0-2D yields -9.9 kcal/mol, and PBE- 
ooD improves the estimate to -11.8 kcal/mol. This agrees 



with our previous findings using an interatomic many- 
body dispersion method [18j . and demonstrates that the 
many-body dispersion energy becomes significant as the 
molecule size increases. 



VI. CONCLUSIONS 

The widely used DFT-D and DFT+vdW methods that 
compute the dispersion energy as a pairwise sum over 
atoms have been derived from the quantum mechanical 
adiabatic connection formula. This derivation puts in- 
teratomic dispersion methods on a firm theoretical basis. 
We have shown that the damping of the dispersion energy 
at short interatomic distances is connected to the spatial 
spread of the dipole moment. We have also demonstrated 
that the non-additive many-body effects beyond the pair- 
wise approximation play an important role for the bind- 
ing energies of dispersion-bound complexes. Moreover, 
given the equivalence between the exact and RPA treat- 
ment of the coupled QHO model, the full many-body dis- 
persion energy can be efficiently computed with a single 
matrix diagonalization. 

There are many avenues remaining for future work, in- 
cluding, for example, (i) different approximations to the 
ACFD formula, (ii) the role of input polarizabilities into 
the ACFD formula, (iii) the role of anisotropy and lo- 
calization in the input polarizabilities, (iv) the role of 
higher multipole moments in the response function, and 
(v) improving the coupling between DFT and the long- 
range dispersion energy. The adiabatic connection for- 
mula provides a powerful framework for the development 
of accurate and efficient approaches for computing the 
correlation energy in general and the dispersion energy 
in particular. 

ACKNOWLEDGMENTS 

All authors acknowledge the European Research Coun- 
cil (ERC Starting Grant VDW-CMAT), and insightful dis- 
cussions with Xinguo Ren and John F. Dobson. 



[1] K. E. Riley, M. Pitonak, P. Jurecka, and P. Hobza, Chem. 

Rev. 110, 5023 (2010) 
[2] F. O. Kannemann and A. D. Becke, J. Chem. Theory 

Comput. 6, 1081 (2010) 
[3] S. Grimme, Comput. Mol. Sci. 1, 211 (2011) 
[4] V. R. Cooper, L. Kong, and D. C. Langreth, Physics 

Procedia 3, 1417 (2010) 
[5] S. Steinmann and C. Corminboeuf, J. Chem. Theory 

Comput. 7, 3567 (2011) 
[6] A. Tkatchenko et al, MRS Bull. 35, 435 (2010) 
[7] G. Mercurio, E. R. McNellis, I. Martin, S. Hagen, 

F. Leyssner, S. Soubatch, J. Meyer, M. Wolf, P. Tegeder, 

F. S. Tautz, and K. Reuter, Phys. Rev. Lett. 104, 036102 



(2010) 

[8] N. Atodiresei, V. Caciuc, P. Lazic, and S. Bfiigel, Phys. 

Rev. Lett. 102, 136809 (2009) 
[9] K. Tonigold and A. Gross, J. Chem. Phys. 132, 224701 

(2010) 

[10] D. Stradi, S. Barja, C. Diaz, M. Garnica, B. Borca, J. J. 

Hinarejos, D. Sanchez-Portal, M. Alcami, A. Arnau, A. 

L. Vazquez de Parga, R. Miranda, and F. Martin, Phys. 

Rev. Lett. 106, 186102 (2011) 
[11] M. Dion, H. Rydberg, E. Schroder, D. C. Langreth, and 

B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004) 
[12] K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and 

D. C. Langreth, Phys. Rev. B 92, 081101 (2010) 



9 



O. A. Vydrov and T. Van Voorhis, J. Chem. Theory 
Comput. 8, 1929 (2012) 

A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 
073005 (2009) 

N. Marom, A. Tkatchenko, M. Rossi, V. V. Gobre, 
O. Hod, M. Scheffler, and L. Kronik, J. Chem. Theory 
Comput. 7, 3944 (2011) 

S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. 
Phys. 132, 154104 (2010) 

O. A. von Lilienfeld and A. Tkatchenko, J. Chem. Phys. 
132, 234109 (2010) 

A. Tkatchenko, R. A. DiStasio, Jr., R. Car, and M. Schef- 
fler, Phys. Rev. Lett. 108, 236402 (2012) 
R. A. DiStasio Jr., O. A. von Lilienfeld, and 

A. Tkatchenko, Proc. Natl. Acad. Sci. USA 109, 14791 
(2012) 

F. London, Z. Physik. Chem. (Leipzig) Bll, 222 (1930) 
H. Margenau, Rev. Mod. Phys. 11, 1 (1939) 
O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 
4274 (1976) 

D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 
(1977) 

S. L. Adler, Phys. Rev. 126, 413 (1962) 

N. Wiser, Phys. Rev. 129, 62 (1963) 

D. Bohm and D. Pines, Phys. Rev. 92, 609 (1953) 

M. Cell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 

(1957) 

M. Dion, H. Rydberd, E. Schroder, D. C. Langreth, and 

B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004) 

M. Fuchs and X. Gonze, Phys. Rev. B 65, 235109 (2002) 
F. Furche and T. van Voorhis, J. Chem. Phys. 122, 
164106 (2005) 

F. Furche, J. Chem. Phys. 129, 114105 (2008) 

G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. 
Chem. Phys. 129, 231101 (2008) 

B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. 
Chem. Phys. 130, 081105 (2009) 

J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. 

Angyan, Phys. Rev. Lett. 102, 096404 (2009) 

J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008) 

J. Harl and G. Kresse, Phys. Rev. Lett. 103, 056401 

(2009) 

D. Lu, Y. Li, D. Rocca, and G. Galli, Phys. Rev. Lett. 
102, 206411 (2009) 

J. F. Dobson and J. Wang, Phys. Rev. Lett 82, 2123 
(1999) 

M. Rohlfing and T. Bredow, Phys. Rev. Lett 101, 266106 
(2008) 

X. Ren, P. Rinke, and M. Scheffler, Phys. Rev. B 80, 
045402 (2009) 

L. Schimka, J. Harl, A. Stroppa, A. Griineis, M. Mars- 
man, F. Mittendorfer, and G. Kresse, Nature Materials 
9, 741 (2010) 



[42] Y. Li, D. Lu, H.-V. Nguyen, and G. Galli, J. Phys. Chem. 

A 114, 1944 (2010) 
[43] X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. 

Rev. Lett. 106, 153003 (2011) 
[44] D. Lu, H.-V. Nguyen, and G. Galli, J. Chem. Phys. 133, 

154110 (2010) 

[45] X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, 

A. Sanfilippo, K. Router, and M. Scheffler, New. J. Phys 

14, 053020 (2012) 
[46] J. F. Dobson, in Topics in Condensed Matter Physics, 

edited by M. P. Das (Nova (New York), 1994) p. 121 
[47] A. Krishtal, K. Vannomeslaeghe, D. Geldof, C. V. 

Alsenoy, and P. Geerlings, Phys. Rev. A 83, 024501 

(2011) 

[48] A. Tkatchenko, D. Alfe, and K. S. Kim, J. Chem. Theory 
Comput. (2012) 

[49] B. M. Axilrod and E. Teller, J. Chem. Phys. 11, 299 
(1943) 

[50] A. J. Stone, The theory of intermolecular forces (Oxford 

University Press, 1996) 
[51] A. Mayer, Phys. Rev. B 75, 045407 (2007) 
[52] A. Koide, J. Phys. B: Atom. Molec. Phys. 9, 3173 (1976) 
[53] S. Grimme, S. Ehrlich, and L. Goerigk, J. Comput. 

Chem. 32, 1456 (2011) 
[54] W. L. Bade, J. Chem. Phys. 27, 1280 (1957) 
[55] A. G. Donchev, J. Chem. Phys. 125, 074713 (2006) 
[56] M. W. Cole, D. Velegol, H.-Y. Kim, and A. A. Lucas, 

Mol. Simul. 35, 849 (2009) 
[57] J. Cao and B. J. Berne, J. Chem. Phys. 97, 8628 (1992) 
[58] H. Berthoumieux and A. C. Maggs, Europhys. Lett. 91, 

56006 (2010) 

[59] X. Ren, P. Rinke, C. Joas, and M. Scheffler, J Mater. Sci. 
47, 7447 (2012) 

[60] J. Applequist, K. R. Sundberg, M. L. Olson, and 
L. Weiss, J. Chem. Phys. 70, 1240 (1979) 

[61] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996) 

[62] F. L. Hirshfeld, Theor. Chim. Acta 44, 129 (1977) 

[63] T. Takatani, E. G. Hohenstein, M. Malagoli, M. S. Mar- 
shall, and C. D. Sherrill, J. Chem. Phys. 132, 144104 
(2010) 

[64] V. Blum et ai, Comp. Phys. Commun. 180, 2175 (2009) 
[65] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. 

Phys. 105, 9982 (1996) 
[66] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 

(1999) 

[67] A. D. Becke, J. Chem. Phys. 88, 2547 (1988) 
[68] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 
(1988) 

[69] P. Jurecka, J. Sponer, J. Cerny, and P. Hobza, Phys. 
Chem. Chem. Phys. 8, 1985 (2006) 



